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ABSTRACT 

If gravitational clustering is a hierarchical process, the present large-scale structure of 
the galaxy distribution implies that structures on smaller scales must have formed at high 
redshift. We simulate the formation of small-scale structure (average cell mass: Amj = 
1O 4 ' 2 M0) and the evolution of photoionized gas, in the specific case of a CDM model with 
a cosmological constant. The photoionized gas has a natural minimal scale of collapse, 
the Jeans scale (m& ; j ~ 1O 9 M ). We find that low column density (N HI < 10 14 cm" 2 ) 
lines originate in regions resembling Zel'dovich pancakes, where gas with overdensities in 
the range 3 — 30 is enclosed by two shocks but is typically re-expanding at approximately 
the Hubble velocity. However, higher column density (Nhi > 10 15 cm -2 ) lines stem from 
more overdense regions where the shocked gas is cooling. We show that this model can 
probably account for the observed number of lines, their distribution in column density 
and b-parameters, as well as the cloud physical sizes as observed in gravitationally lensed 
quasars. We find a redshift evolution that is too steep; however, this may be due to 
insufficient dynamical range in the simulation or because the specific model (CDM+A) is 
incorrect. The model predicts that high signal-to-noise observations should find systematic 
deviations from Voigt profiles, mainly in the form of broad wings in the line profiles, 
and that a fluctuating Gunn-Peterson effect will be detected, which can be modeled as a 
superposition of weak lines with a wide range of b-parameters. 

Cosmology: early universe -Cosmology: large-scale structure of Universe - hydrodynamics 
- quasars: absorption lines 
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1. INTRODUCTION 

The physical origin of the Ly a forest in the QSO absorption spectra has remained 
a mystery. The essential, difficult fact is that the clouds have "size" ~ 10 — 30kpc and 
are typically separated by ~ 2,000kpc along the line of sight, but produce most of the 
absorption. Thus, the absorption per unit length within a cloud must exceed the absorption 
per unit length in the intercloud medium by a factor of more than 10 2 . Because the neutral 
hydrogen is proportional to p 2 , the clouds must then be overdense by a factor > 10. This 
has led to extensive discussions of the possible confinement mechanism, with gravity and 
external gas pressure being proposed, as well as free expansion. We follow up the scenario 
of gravitational collapse on subgalactic scales for the origin of the Ly a clouds (Ikeuchi 1986; 
Rees 1986; Bond, Szalay, & Silk 1988). We go beyond the previous qualitative models, by 
performing direct, high resolution, numerical hydro simulations of gravitational collapse 
on small scales. We find that we do, in fact, reproduce most of the observed properties 
of the Ly a forest in a simulation of the CDM+A model, although, as we shall see, the 
physical picture is different from that anticipated in prior semi-analytic calculations. 

If this theory of the Ly a forest is correct, then the observations of Ly a extinction in 
quasars will provide a new, direct probe to the initial density fluctuations in the universe. 
Such observations have three advantages when compared to local (z = 0) observations of 
galaxies for the assessment of cosmological models: (1) they are unbiased in that they 
probe random lines of sight; (2) they can reach to relatively high redshift (z > 4), where 
the contrasts between competing models are greatest; (3) they provide a large database 
with ~ 10 2 lines/quasar xlO 3 quasars ~ 10 5 lines as the potential data base. 

2. NUMERICAL MODELING 

We use the new shock-capturing cosmological hydrodynamic code described in Ryu 
et al. (1993) called Total Variation Diminishing (TVD) with addition of atomic processes 
within a primeval plasma of (H,He). We calculate self-consistently the average background 
photoionizing radiation field as a function of frequency. Ionization equilibrium of the 
abundances of HI, HII, Hel, Hell, and Helll is required at every timestep, but the time 
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evolution of the temperature in every cell is solved explicitly, using the same heating and 
cooling terms as in Cen (1992). 

We model galaxy formation as in Cen & Ostriker (1992, 1993a, b). The material 
turning into galaxies is assumed to emit ionizing radiation, with two types of spectra: 
one characteristic of star formation regions and the other characteristic of quasars, with 
efficiencies (i.e., the fraction of rest-mass energy converted into radiation) of ejjv,* = 
5 X 10 -6 , and ejjv,Q = 6 X 10 -6 , respectively. The derived intensity of the radiation field 
at the Lyman limit is J = (2.3,3.4,1.8) X 10~ 21 ergcm~ 2 sec~ 1 hz~ 1 sr~ 1 at z = (2,3,4), 
slightly higher than suggested by observations of the proximity effect (Murdoch et al. 
1986; Bajtlik, Duncan, & Ostriker 1988; Lu et al. 1991). 

We adopt a CDM+A model with the following parameters: O = 0.4, = 0.0355 
(cf. Walker et al. 1991), A = 0.6, h = 0.65 and <j 8 = 0.79 (COBE normalization; Kofman, 
Gnedin, & Bahcall 1993). This model is consistent with all observational tests that we 
know of. Nevertheless, we have no reason to believe that it is, in fact, correct. Our box size 
is 3/i _1 Mpc with N = 288 3 cells and 144 3 dark matter particles. The cell size is 10.4/i _1 kpc 
(comoving) corresponding to a baryonic mass of 1.7 X 1O 4 M0, with true spatial resolution 
slightly worse than this. At redshift three, the Jeans length, Aj = (vrc 2 / 'Gptot) 1 ^ 2 for 
c s = v rms = lOkrns -1 , is equal to 400/i _1 kpc in comoving units, or 40 cells, so the smallest 
scales of initial collapse for the photoionized g well resolved. 

3. PHYSICAL PROPERTIES OF THE COLLAPSED STRUCTURES 

Figure 1 shows contour plots of four quantities projected accross a slice of the sim- 
ulation of thickness 0.5/i _1 Mpc: total gas density neutral hydrogen column density 
accross the slice (p_fjj), dark matter density (pd) } and gas temperature (T). We see that 
high column density Ly a clouds (iVjjj > 10 14 cm -2 ) are typically isolated regions, while 
lower column density clouds are generally in connected filaments and sheets. As expected, 
photoionization equilibrium prevails so that T ~ 10 4 K and njji ot n 2 Htot . 

Next, we take a random row along the simulation, and examine the physical structures 
that are intercepted. The upper panel in Figure 2 shows T (solid line) and p& (thick dotted 
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Figure 1. Coutour plots of four quantities in a slice of 1.5 X 1.5 X 0.5/i~ 3 Mpc 3 ( comoving) at redshift 
three: density, neutral hydrogen column density, dark matter density and baryonic temperature. Contour 
levels are 10 ' 1 ' 2 for the total gas density and dark matter density, each in terms of its own global mean, 
and lO^^cm 2 for the neutral hydrogen column density (z = 0,1,2,...). The contour levels for 
temperature are 

line) along such a row, as a function of the spatial coordinate, given as v = x H(H = 
512 h kms _1 Mpc _1 at z = 3). The thin dotted line is the gas pressure (p), in arbitrary 



units. The highest density peak seen resembles the structure of a Zel'dovich pancake 
(Zel'dovich 1970; Sunyaev & Zel'dovich 1972), with two shocks propagating outwards, and 
a minimum of the temperature at the center, where the gas has cooled after being shocked. 
The central minimum exists even in an adiabatic calculation as the shock velocity (and 
thus the post shock temperature) grows with time. 

We can now re-examine the question of the dynamical mechanism that causes the gas 
overdensity, and therefore the absorption lines. The possible mechanisms are: (i) Pressure 
confinement, as proposed by Sargent et al. (1980) and Ostriker & Ikeuchi (1983) does not 
play a role, since the temperature is higher in the density peaks than in the surrounding 
medium, (ii) The inertia of the gas; the gas must then be in free expansion (Ikeuchi & 
Ostriker 1986; Bond, Szalay, & Silk 1988). (iii) Gravity (Rees 1986; Ikeuchi 1986; Ikeuchi 
& Ostriker 1986). (iv) Ram pressure by infalling gas, which has not been proposed so far. 
The latter two processes seem to be dominant. 

To see this we plot in Figure 2b the gravitational acceleration divided by the Hubble 
constant, —jjjf^, as the dashed line, and the total acceleration, — (77 ), as the solid 

line. The dotted line is the peculiar velocity (v p ) along the row; in the linear regime, and 
for = 1, this is given by v p = — In the highest density peak in the figure, the spikes 
in the total acceleration and the discontinuities in the velocity correspond to the two shock 
fronts. Between the shocks, there is a large gravitational acceleration, and the solid line 
shows that this is approximately balanced by the pressure gradient. The ram pressure of 
the infalling gas at the shock is, in this case, about one third of the maximum pressure in 
the system, so the ram pressure and gravitational binding forces are comparable. We have 
examined a total of 30 random rows and find that this situation is typical (although the 
balance between gravitational and pressure forces between shocks is not generally perfect), 
and the gas is typically expanding in comoving coordinates. In other rows, where the 
overdensity is higher than 300, in at least one cell (and Nhi > 10 15 cm -2 ), the gas is 
typically contracting and cooling due to collisional excitation processes. 

We then conclude that, although the Ly a clouds that we simulate arise from gravita- 
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Figure 2. Panel (a) shows the gas temperature (solid line) and gas density (thick dotted line) along such 
a row, as a function of the spatial coordinate, given as V = X H . The thin dotted line is p, in arbitrary 
units. Panel (b) shows the gravitational acceleration divided by the Hubble constant, — J{ ~Jx ' as the 

dashed line, and the total acceleration — tt(~^~ ~\~ as the solid line. The dotted line is the peculiar 

' H ^ ax p ax / ' 1 

velocity along the row. Panel (c) shows the flux reduction from unity for the row featured in Figure 2a 
plotted as the solid line. The spectrum convolved with a gaussian of width = \2,kT(x~) / Uljj )] 1/ ' 2 is 

shown as the dashed curve. The dotted curve is obtained by convolving the dashed curve with a gaussian 
of FW H M = 10 kms 1 to simulate the instrumental resolution, with added gaussian noise having 
S/N = 15 on a sampling bin of 3 kms (e.g., Carswell et al. 1991). 
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tional collapse, as in Bond, Silk & Szalay (1988) and the minihalo model of Rees (1986) and 
Ikeuchi (1986), they differ from this model in three important respects: first, most of the 
absorption lines do not correspond to virialized regions with the gas in hydrostatic equi- 
librium, but rather arise from the infalling gas forming the small-scale structure (Miralda- 
Escude & Rees 1993). Second, the gas of low column density clouds (Nhi < 10 14 cm -2 ) has 
been shocked, and the ram pressure at the shocks plays an important role in maintaining 
the overdense region that causes the absorption line, with the post shock gas expanding as 
in a Zel'dovich pancake. Third, the gas does not need to be in perfect thermal equilibrium 
with the ionizing background, because it is shock-heated and, depending on the intensity 
and spectrum of the ionizing background, the cooling time at the typical densities of the 
clouds can be comparable to the Hubble time (Efstathiou 1992). Shocks will be quasi- 
isothermal with the consequence that large density jump (p' j p) % (7 + l)/(7 — 1) occurs 
in the transition. 

We now use (p/jj, v p ) to calculate the optical depth due to Ly a scattering, if there 
were no thermal motions of the hydrogen atoms: t(V) = tqp {1 + ^/[^(V")]} ^y. 
Here, V is the component of the total velocity of the gas along the line of sight, with 
respect to an arbitrary origin; xi(V) are the set of all the positions along the line of sight 
where the gas moves with velocity V; and Shi is the overdensity of neutral hydrogen. In 
the linear regime, there is only one spatial position x\ for every velocity V", but when 
structures start collapsing velocity caustics appear, and within those x(V) is multivalued. 
The normalization tqp (Gunn-Peterson optical depth) is the optical depth if the medium 
had constant density. Thus, tqp = 4.6 X 10 5 hy (1 + z) 3 / H(z), where is the baryon 
density divided by the critical density, y is the neutral fraction of hydrogen, and H(z) is 
the Hubble constant at redshift z. For the present simulation, tgp = 0.052 at z = 3. 

The flux reduction from unity (e _T ) for the row featured in Figure 2a is plotted as 
the solid line in Figure 2c. The resultant spectrum of convolving the optical depth with a 
gaussian of width b ca [ c (x) = \2kT(x) / run )] 1 ^ 2 is plotted as a dashed curve in Figure 2c. In 
general, the absorption lines arise from density peaks, producing narrow and well-defined 
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features. Velocity caustics (see McGill 1990) are not important in generating lines due to 
thermal broadening. 

We now convolve the dashed curve with a gaussian of FWHM= 10 kms" 1 to simulate 
the instrumental resolution, and add gaussian noise with S/N = 15 (adequate for current 
observations, e.g., Carswell et al. 1991) on a sampling bin of 3 kms -1 , obtaining the 
dotted curve in Figure 2c. We repeat the same procedure on 30 randomly selected rows, 
and detect absorption features and fit Voigt profiles (cf. Webb 1987 and Carswell et al. 
1987) to the simulated spectra in the same way as it is usually done for real data (e.g., 
Rauch et al. 1992). The statistical significance of any spectral features is determined 
using the procedure described by Young et al. (1979), with lines being considered real 
detections when their equivalent width has a probability less than 10 -5 to arise by random 
fluctuation. Whenever \ 2 nas l ess than 1% probability of occurring by chance, additional 
components are included to fit the features as blends. 

The b-parameters of the lines are in the range 15 - SOkms" 1 (b calc = 20kms" 1 ), 
slightly lower than observed, the later being in the range 15 — 60 kms -1 (b i, s ~ 35kms _1 ) 
(cf. Press & Rybicki 1993, Rauch et al. 1992 for values and references to earlier literature). 
There are two effects that could increase the b-parameters: first, the temperature of the 
intergalactic medium (8000K in this simulation) could be higher depending on the epoch of 
reionization and the spectrum of the ionizing background (Miralda-Escude & Rees 1994). 
Second, our small box simulation does not include the large-scale power, which would 
cause broader lines. In addition, there is evidence that the b-parameters are lower at high 
redshift (Williger et al. 1994), in better agreement with our derived values. In any case, 
bulk motions and shock-heating are inherent to our model, and imply that b-parameters 
should be larger than for photoionized gas in thermal equilibrium, as is observed. 

Next, we take a random row along the simulation, and examine the physical structures 
that are intercepted. Next, consider the column density distribution of the clouds. We 
find that the lines that are identified as absorption systems generally correspond well with 
overdense regions in real space. Therefore we can obtain the column density distribution 
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Figure 3. Thin solid and dotted lines show the column density distribution at redshift Z = (3, 2). Solid 
dots indicate the observed number of absorption lines with Njjj > 10 14 cm 2 (Carswell et al. 1991), 
and N H I > 10 13 - 75 cm~ 2 and N H I > 10 14 - 27 cm~ 2 (Rauch et al. 1992). Short dashed lines show 
the range of slopes of the column density distribution around N HI = 10 14 cm- 2 obt ained by different 
observational methods (Rauch et al. 1992; Petitjean et al. 1993; Press & Rybicki 1993). Solid and dotted 
lines are extraplated column density distribution at Z = (3, 2) for J ~ 10 21 . In addition, we show, as 
the thick dashed curve, the result from a 144 3 simulation of the same model scaled to J = 10 21 , to 
show the numerical resolution effect. 

approximately by taking, through each column of the simulation at z = 3 (a total of 
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288 X 288 X 3), all the connected regions of cells with an overdensity larger than unity, 
considering each such region as a "cloud". The result is shown in Figure 3, at redshifts 
z = (3,2), as the thin solid and dotted lines respectively. Also shown as a solid square 
is the observed number of absorption lines with Nhi > 10 17 ' 3 cm -2 from Sargent et al. 
(1989), and solid dots from Rauch et al. (1992), and an open circle from Petitjean et al. 
(1993). The short dashed lines show the range of slopes of the column density distribution 
around Nhi = 10 14 cm -2 obtained by different observational methods (Rauch et al. 1992; 
Petitjean et al. 1993; Press & Rybicki 1993). The column density distribution we obtain is 
close to but somewhat steeper than observed. The number of lines with N > 10 13 ' 75 cm -2 
per unit redshift is smaller than observed by a factor of ~ 2 (Rauch et al. 1992) 

Several effects could alter our derived column density distribution. First, if the inten- 
sity of the ionizing background is lower than we have computed, the curves in Figure 3 are 
all shifted horizontally to higher column densities. For example, if J = 10 -21 , in agree- 
ment with the observed proximity effect, then the number of lines at Nhi > 10 14 cm -2 
is increased to the observed value [shown as thick solid and dotted lines for redshift 
z=(3,2)] and the slope of the distribution shallower at the same column density (around 
Nhi = 10 14 cm -2 ), in better accord with the observations. Second, due to the limited 
resolution of the simulation, the number of systems with high column densities is probably 
underestimated (illustrated by the thick dashed curve from a 144 3 simulation of the same 
model scaled to J = 10 -21 ). Finally, due to the small size of our box, large-scale power 
is missed, and systems arising from structures on large scales are not taken into account. 
Correcting this should produce more systems with higher column densities, decreasing the 
slope of the distribution at Nhi = 10 14 cm -2 . 

We now examine the redshift evolution of the clouds. If the redshift evolution is 
parameterized as N(z) oc (1 + z) 7 , we find, for lines with Nhi > 10 14 cm -2 , 7 = (6.9,4.1) 
over the ranges z = (4 — > 3, 3 — > 2), respectively, and 7 = (4.9, 1.8) over the same ranges 
for Nhi > 10 16 cm -2 . The evolution of the number of lines with redshift is sensitive to 
the variation of J. If J is constant from redshift four to two, the same values of 7 change 
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to (3.5,4.5) and (2.5,2.2), respectively. The redshift evolution appears to be steeper than 
observed (e.g., Lu et al. 1991), although the fact that the evolution is less pronounced at 
larger column densities agrees with observations. The rapid destruction of the clouds is 
presumably due to their expansion and to merging into larger structures. Also, different 
cosmological models will predict different rates of evolution. 

We can also see from Figure 1, from the contours of the neutral column density, that 
lines with Nhi ~ 10 14 cm -2 have transverse sizes of ~ 20/i _1 kpc, which agrees with the 
observational constraints on the physical size of the clouds from gravitationally lensed 
quasars (Foltz et al. 1984; Smette et al. 1992). Thus, we see that the structure found in 
the simulation does in fact roughly correspond to the well studied "Lyman alpha forest". 

Our model can make predictions for the high signal-to-noise observations of the Ly a 
forest with the new generation of large telescopes. We have analyzed the same 30 random 
spectra assuming S/N = 60 on a sampling bin of 3 kms -1 , and resolution of 7 kms -1 . Of 
60 features detected, 46 are normal lines. Seven show broad wings which are modeled as 
blends of two systems at almost the same wavelength, but very different b-parameters; this 
is caused by high- velocity gas infalling into the system, and a distribution of temperatures 
(the strongest line in Fig. 2c is an example). Another 7 are "Gunn-Peterson features", 
with Nhi < 10 12 ' 7 cm -2 and b > 30 kms -1 , caused by small fluctuations in the IGM. When 
the signal-to-noise is reduced to 15, only two lines with broad wings could be detected, 
and no Gunn-Peterson feature was found, out of 37 detected lines. 

Deviations from Voigt profiles are inherent to our model. When attempting to fit 
them with blends, they should cause a large peak in the spatial correlation of the lines at 
distance of order ~ 6, which increases with signal-to-noise. A second prediction is that a 
fluctuating Gunn-Peterson effect is present between the lines, but this can be modeled as 
low column density lines. 

5. CONCLUSIONS 

We have shown that the growth of structure on subgalactic scales by gravitational 
collapse can probably account for the observed Ly a forest. The low column density ab- 
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sorption lines are produced by photoionized gas in regions resembling Zel'dovich pancakes, 
where gas is shock-heated and is maintained at an overdensity of 3 to 30 by a combination 
of gravity and ram pressure from infalling gas. The geometry of the low column density 
"clouds" varies from sheet-like to filamentary to quasi-spherical. 

Ly a clouds should therefore be considered as probes to the formation of structure 
and to the primordial density fluctuations of the universe, in addition to galaxies. The 
calculation of the density and temperature distribution of the neutral gas that we have 
done in our cosmological simulation is still subject to several uncertainties, due to the 
limited resolution of the simulation, the exclusion of large-scale power, and the poorly 
known ionization history of the intergalactic medium and the spectrum of the ionizing 
radiation. Some of these problems will be improved by future simulations with a higher 
dynamic range. Given the wealth of observational data that is available on the quasar 
absorption systems, we may expect that their study will provide sensitive new constraints 
on theories of structure formation. 

This work is supported in part by NASA grant NAGW-2448 and NSF grants AST91- 
08103 and ASC-9318185, and the W. M. Keck Foundation. We thank Martin Rees for 
stimulating discussions, and Bob Carswell and John Webb for providing us with their 
profile fitting software. The simulation was performed on Convex-3880 in NCSA with time 
allocation provided by M. Norman. 
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